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ABSTRACT 

We utilize magnetohydro dynamic (MHD) simulations to develop a numerical model for GMC-GMC 
collisions between nearly magnetically critical clouds. The goal is to determine if, and under what 
circumstances, cloud collisions can cause pre-existing magnetically subcritical clumps to become su¬ 
percritical and undergo gravitational collapse. We first develop and implement new photodissociation 
region (PDR) based heating and cooling functions that span the atomic to molecular transition, cre¬ 
ating a multiphase ISM and allowing modeling of non-equilibrium temperature structures. Then in 
2D and with ideal MHD, we explore a wide parameter space of magnetic field strength, magnetic field 
geometry, collision velocity, and impact parameter, and compare isolated versus colliding clouds. We 
find factors of ~ 2 — 3 increase in mean clump density from typical collisions, with strong dependence 
on collision velocity and magnetic field strength, but ultimately limited by flux-freezing in 2D geome¬ 
tries. For geometries enabling flow along magnetic field lines, greater degrees of collapse are seen. 
We discuss observational diagnostics of cloud collisions, focussing on 13 CO(J=2-l), 13 CO(J=3-2), 
and 12 CO(J=8-7) integrated intensity maps and spectra, which we synthesize from our simulation 
outputs. We find the ratio of J=8-7 to lower- J emission is a powerful diagnostic probe of GMC 
collisions. 

Subject headings: ISM: molecular clouds — ISM: magnetic fields — methods: numerical 


1. INTRODUCTION 

Understanding star formation is a key astrophysical 
problem, with many fundamental questions still unre¬ 
solved. In particular, what mechanisms drive or inhibit 
the star formation process? Studying the evolution of gi¬ 
ant molecular clouds (GMCs) within the diffuse interstel¬ 
lar medium and the formation of prestellar clumps and 
cores within GMCs is complicated because of the large 
range of densities (nn ^ 1 cm -3 to ~ 10 6 cm -3 ), length 
scales (^kpc to ~0.1 pc), and timescales 10 8 yr for 
galactic orbits to ~ 10 5 yr for core dynamical timescales) 
involved, as well as the nonlinear effects of self-gravity, 
thermal, magnetic, and turbulent pressures, large-scale 
motions such as galactic shear or collisional converging 
flows, radiation, chemistry, and feedback. Additionally, 
the initial conditions are uncertain and boundary con¬ 
ditions are poorly constrained. The final conditions, 
gleaned through observations, suggest that star forma¬ 
tion is highly clustered and localized, with relatively high 
local effi cien c y with in clusters (e.g., Lada & Lada 2003 


Gutermuth et al.]|2009 ). However, overall star formation 


is slow and inefficient, with only a few percent of total 
gas forming stars over local dynamical timescales (e.j 


Zuckerman & Evans||1974| IKrumholz & Tan||2007| |Da| 
Rio et al]2014| ). -"- M -"- 


W ith regard to the importance o f magnetic fields (see 


Crutcher||2012 Li et al.||2014 ), 


there have been two 


e -g- 

main views. Strong-field models propose relatively long 
GMC lifetimes in which magnetic fields play important 
roles in controlling formation and evolution of the clouds. 
In these models, non-star-forming clouds are initially 
subcritical, i.e., their magnetic fields are strong enough 
to prevent gravitational collapse. Weak-field models have 
GMCs as intermittent phenomena with short lifetimes 
10 6 yr) and turbulent flows controlling the formation 
of clouds, clumps and cores. These models posit mag¬ 
netically supercritical masses, i.e., the magnetic pressure 
alone is too weak to support against gravity. 

Zeeman measurements show that mass-to-magnetic 
flux ratios (M/4>) are approxim ately critical to slightly 
supercritical in molecular clouds ([Crutcher 1999||Troland 
Crutcher|2008||Crutcher|2012||Li et al.|2(J14|). If GlVICs 
are partially stabilized by magnetic fields, then this may 
increase their lifetimes to > 20 Myr timescales, which 
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are then comparable to GMC-GMC collision ti mes, es¬ 
pecially for clouds inside the solar circle (Gammie et al 


1991 Tan 


1|2QQQ| [Tasker fe Tan||2009| [Dobbs etaU|2015[). 

mservational evidence tor frequent GlVlC col- 


Indirect 

lisions comes from the near random orientations of ' 


jected angular momentum vectors of GMCs (Rosolowsky 

et al. 2003 

Koda et al. 2006; Imara & Blitz||2011; Imara 

eta! 2011 

). 


anism for injecting turb ulent energy into GMCs (Tan 


2000 Tan et al. 


2013), with the energy being ex¬ 


tracted from galactic orbital motion. Other mechanisms 
for injecting tu r bulence involve star for mation feedback 
(|Matzner||2007| |Goldbaum et al.||2011|). Without such 
replenishment, turbule nce is expected to decay within 
about a cr ossing time (Mac Low et al. 1998 Ostriker 
et al. 1999). 


By producing dense gas, compressed in shocks, GMC- 
GMC collisions may bejm important trigger of star clus¬ 
ter formation (jScoville et al.|1986). If the majority of star 
formation is initiated by this process, then a model of 
shear-mediated cloud collisions can naturally exp lain the 
obser ved dynamical Kennicutt-Schmidt relation (Kenni- 


cutt 1998) in which a roughly constant fraction of gas, 
e or b — 0.04 is converted into stars every local orbital 
time ([Tan 2000, 2010 Suwannajak et al. 2014). Note 
that this mechanism of creating star-forming molecular 
clumps from localized compressed regions of pre-existing 
GMCs, is different from that proposed for creating molec- 



Heitsch et al. 20061 

van Loo et al. 2007 

et al. 

2009; 

jjvan Loo et al. 

2010p 


Heitsch 


number of previous studies. Habe & Ohta (1992) per¬ 
formed 2D axisymmetric SPH simulations of head-on col¬ 
lisions of non-identical clouds. These collisions produced 
a bow shock which disrupted the larger cloud while com¬ 
pressing the smaller cloud. This compression could lead 
to gravitational instability for the smaller cloud even if 
its initial mass was below the Jeans mass. 


Klein & Woods (1998) presented 2D AMR hydrody¬ 
namics simulations of homogeneous cloud collisions. The 
collisions resulted in bending mode instabilities creating 
large aspect ratio filaments. With surface perturbations, 
the merged cloud system became highly asymmetrical 
and highly inhomogeneous with islands of high density 
surrounded by low density regions. 


Anathpindika (2009) performed a series of 3D SPH 
simulations which investigated the gravitational stability 
of post-shock compressed slabs resulting from molecular 
cloud collisions. Additionally, sheared collisions result in 
non linear thin shell instabilities and Kelvin-Helmholtz 
instabilities. 

More recently, 


Takahira et al. 


(2014) performed 3D 


hydrodynamic simulations with AMR showing core for 
mation occurring from a GMC collision interface. They 
found that faster collision velocities formed a greater 
number of cores, but core growth was predominantly via 
accretion in the shock front, with slower shocks being 

favored for making larger cores. _ 

In terms of MHD collision studies, Kortgen fc Baner^l 


jee (2015) investigated molecular cloud formation and 
the transition from magnetically sub- to supercritical HI 
clouds via converging magnetized flows. Even with mag¬ 


netic diffusion effects, they found that cylindrical flows 
created no magnetically supercritical regions and star 
formation is strongly suppressed for even relatively low 
initial magnetic field strengths. 

This paper, the first of a series, explores the process of 
magnetized cloud-cloud collisions and its effect on indi¬ 
vidual GMC and clump scales. Here we restrict analysis 
to a parameter space exploration with 2D simulations 
of simplified cloud geometries, including an embedded 
clump: formally colliding infinite cylinders, which can 
approximate collisions of spheroidal clouds. 

^explains the fiducial set-up and various simulation 
and analytic methods employed. m describes new heat¬ 
ing and cooling functions that we have developed for 
this project. ^4] describes the set up and subsequent 
results of exploring the following parameters: magnetic 
field strength, magnetic field orientation, collision veloc¬ 
ity, and impact parameter. discusses predictions of 
observational diagnostics of shocks. Discussion and con¬ 
clusions follow in ^6j 


2. NUMERICAL MODEL 
2.1. Initial Conditions 

For our default initial conditions we use typical ob¬ 
served values of Galactic GMC and ISM properties. 
GMCs are conventionally defined as having masses 
> 10 4 M©. They 1 
v - 100 M© pc -2 (e 


mean mass surfaci 

e densities 

McKee & Ostriker 

2007 

Tan 


100 cm“^7 although clumps and cores within the clouds 
have densities that can be orders of magnitude larger. 

GMCs have internal velocity dispersions that are sim¬ 
ilar to virial velocities, typically several km/s, which is 
much larger than the ^0.2 km/s sound speeds of ^ 10 K 
gas. Supersonic turbulence and self-gravity are thought 
to be two important processes that help give rise to the 
hierarchical density structures seen in GMCs. However, 
these structures may also be regulated by magnetic fields. 

The magnetic field of the local diffuse ISM background 
is 6 =b 2 /iG ( Beckp OOl). If a random, uniform distribu¬ 
tion of field strengths is assumed up to a maximum value, 
Smax, Zeeman measurements reveal that this maximum 
magnetic field value measured within molecular clouds, 
clumps and cores with nn > 300 cm -3 sc ales as = 

^o(lRi/300cm -3 ) 0 - 65 , where = 10 /iG (Crutcher et al. 
2010). At lower densities, L> max = Bq = 10 /iG, in- 
dependent of density. We will henceforth refer to this 
as the “Crutcher relation”. Thus, for nu = 10 3 cm -3 , 
L?max — 22 /iG. 

Observed random velocities of Galactic GMCs are ap¬ 
proxi mately 5-7 km s -1 (e.g., |Liszt et al. 1984 Stark 
1984). However, interaction velocities between colliding 
GMCs are likely to be s et by the shear velocity at 1 to 2 
tidal radii of the clouds (|Gammie et al.|1991| |Tan|2000|), 
which can be several times larger. 

Given the 2D nature of the simulations of this paper, 
the modeled structures can be considered as filaments 
extending perpendicular to the simulation domain. We 
follow “clouds”, i.e., “GMCs,” with a uniform density 
of ^h,gmc = 100 cm -3 . Although the clouds are, in 
principle, cylinders of infinite extent, we assume that the 
clouds have a finite mass, i.e., Mqmc = 10 5 M©. A 
mass surface density of E = 100 M© pc -2 integrated 
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TABLE 1 

GMC and Clump Properties 



M to f\ 

( 10 5 Me) 

n H 

(cm- 3 ) 

R 

(pc) 

B cri M 

(/xG) 

Ambient 

- 

10 

- 

10 

GMC f\ 

1.78 

100 

23.8 

42.0 

GMC 2 

0.56 

100 

11.9 

13.9 

Clump 

0.10 

1000 

5.64 

66.0 


a Masses are for an equivalent spherical cloud. 
b Critical B-field strengths are listed for the GMCs and the clump, 
along with the fiducial ambient field strength. 

C GMC 1 includes a clump, but properties listed here are for non¬ 
clump material within the cloud 

along the cloud axis, then, gives a typical cloud radius 
Rgmc = 17.8 pc. We set the radius of the first cloud, 
i.e., Cloud 1, to R\ = 23.8 pc and of the second one, i.e., 
Cloud 2, to R 2 = 0.5i?i = 11.9 pc. 

GMCs are structured, containing dense clumps. When 
a collision occurs, the effect of this collision on pre¬ 
existing clumps may be the most important for triggering 
star formation. We therefore introduce an idealized em¬ 
bedded, overdense clump into Cloud 1. The clump has 
a uniform density of tzh, c i = 1000 cm 3 , i.e., 10 x over- 
dense compared to the GMC, and a radius of 5.6 pc. We 
position the clump off center at (x,y) = (0.5Ri,0) (see 
FigJTJ) . The properties of our clouds and clump are listed 
in Table Q~| 



Fig. 1.— Basic cloud collision setup. GMC 1 (left cloud) has 
radius R \. It includes an embedded clump with radius R c \ located 
at a distance of one half-radius to the edge. GMC 1 collides with 
GMC 2, a uniform cloud with a radius R 2 = Ri/2. The clouds 
are initially separated by a distance that is changed depending on 
relative velocity so that collisions occur at the same time. 


For typical molecular cloud temperatures, ^ 10 — 20 K, 
such GMCs and clumps are not thermally supported 
against gravitational collapse. For example, if the clouds 
are considered as long filaments in the direction perpen¬ 


dicular to the simulation plane, then the mass per unit 
length for GMC 1, mi = 6150 M© pc -1 , far exceeds 
the critical line mass for a cylindrical cloud, given by 
crit = 2 c 2 JG 1 whic h is ^ 20 M© pc -1 for a cold 10 K 


cloud (Ostriker 1964). 

The inclusion oF magnetic fields helps to stabilize the 
clouds. We vary the direction, i.e., parallel, perpendic¬ 
ular and oblique to the cylindrical axis of the cloud, 
and the magnitude of the field. The field strengths are 
detailed in the results section, but are of the order of 


10 — 100 /iG, given observed field strengths (Crutcher 
2012 ). 


Internal cloud turbulence is another mechanism that 
may help support GMCs. To separate out the effects 
of magnetic fields from turbulence, in this paper we do 
not initialize the clouds with any turbulence, deferring 
this to Paper II, which also extends the dimensionality 
to 3D. However, turbulent motions are generated by the 
GMC-GMC collision, which may then provide additional 
support to the clouds. 

The ambient medium in which the clouds are embed¬ 
ded, representative of the atomic cold neutral medium, 
is set to have nn,o = 10 cm -3 , with a magnetic field of 
Bq = 10/iG. The default relative collision velocity of the 
clouds is set to be f re i = lOkms -1 , with variations from 
5 to 40 km s _1 . The default impact parameter, 6, of the 
collision is set to zero, i.e., an on-axis collision, but some 
cases with b = 0.5,1,1.5 iii are also explored. The sur¬ 
rounding medium, which we consider to be a co-moving 
atomic envelope around the GMC, is also colliding. Thus 
in terms of the simulation domain, half the box is moving 
with +v re \/2 and the other half has — v re \/2. However, at 
faster velocities > 20km s -1 we sometimes notice modest 
effects of numerical viscosity on clump properties and so 
also run simulations in the velocity frame of Cloud 1. 

A summary of key parameters in all runs performed in 
this paper is listed in Table |2j Velocities denoted with a 
indicate models run in the frame of Cloud 1. 

2.2. Numerical Code 

The numerical code is a modified version of the Adap- 
tive Mesh Refinement (AMR) code Enzo 2.0 (Bryan & 
Norman|1997||Bryan|1999||0’Shea et al.|2004|). To solve 
the magnetohydro dynamical equations, we use a 2nd- 
order Runge-Kutta temporal update of the conserved 
variables with the Local Lax-Friedrichs (LLF) Riemann 
solver and a piecewise linear reconstruction method. To 
ensure the solenoidal constraint on the magnetic field, the 
divergenc e cleaning algorithm of Dedner et al. (2002) is 
adopted ( Wang fe Abel||2008 ). 


We implemented heating and cooling functions in the 
code that describe both atomic and molecular heating 
and cooling processes (see ^3] for details). These func¬ 
tions take into account a density ve rsus col umn d ensity 


extinction relation similar to that of Van Loo et al. (2013 
henceforth, VLBT2013). As the temperature""of the gas 
needs to be calculated accurately, we use a “dual energy 
formalism” by solving the internal energy equation as 
well as the total energy equation. The temperature is 
then determined from the internal pressure when mag¬ 
netic and kinetic energy together exceed 0.999 the total 
energy, and from the total energy otherwise. 

To track the evolution of properties of the clump, we 
use a scalar value to differentiate between gas outside 
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TABLE 2 

Summary of Simulation Run Parameters 


Run 

™H,0 

™H,1 

™H,cl 

nu ,2 

Bo 

Bi 

B c , 

b 2 

^rel 

b 

0 

;m-3) (. 

cm- 3 ) 

(cm- 3 ) 

(cm- 3 ) 

(/xG) 

(/xG) 

(flG) 

(^G) 

(km/s) 

( Ri ) 

0. Out-of-plane Fields (0,0 ,B Z ) 

0. Isolated Cloud 










O.A.O 

10 

100 

- 

- 

(0,0,10) 

(0,0,27.9) 

~ 

- 

- 

- 

0.A.1 

10 

100 

- 

- 

(0,0,0) 

(0,0,0) 

- 

- 

- 

- 

0.A.2 

10 

100 

- 

- 

(0,0,10) 

(0,0,10) 

- 

- 

- 

- 

0.A.3 

10 

100 

- 

- 

(0,0,10) 

(0,0,20) 

- 

- 

- 

- 

0.A.4 

10 

100 

- 

- 

(0,0,10) 

(0,0,40) 

" 

- 

- 

" 

1. Out-of-plane Fields (0,0, B z ) 

l.B. Isolated Cloud with Clump 








1.B.0 

10 

100 

1000 

- 

(0,0,10) 

(0,0,40) 

(0,0,65) 

- 

- 

- 

l.B.1.1 

10 

100 

1000 

- 

(0,0,10) 

(0,0,10) 

(0,0,65) 

- 

- 

- 

l.B.1.2 

10 

100 

1000 

- 

(0,0,10) 

(0,0,20) 

(0,0,65) 

- 

- 

- 

l.B.1.3 

10 

100 

1000 

- 

(0,0,10) 

(0,0,65) 

(0,0,65) 

- 

- 

- 

l.B.2.1 

10 

100 

1000 

- 

(0,0,10) 

(0,0,40) 

(0,0,40) 

- 

- 

- 

l.B.2.2 

10 

100 

1000 

- 

(0,0,10) 

(0,0,40) 

(0,0,55) 

- 

- 

- 

l.B.2.3 

10 

100 

1000 

- 

(0,0,10) 

(0,0,40) 

(0,0,75) 

- 

- 

- 

l.C. Cloud Collision 










1.C.0 

10 

100 

1000 

100 

(0,0,10) 

(0,0,40) 

(0,0,65) 

(0,0,13.2) 

10 

0 

l.C.1.1 

10 

100 

1000 

100 

(0,0,10) 

(0,0,40) 

(0,0,65) 

(0,0,13.2) 

5 

0 

l.C.1.2 

10 

100 

1000 

100 

(0,0,10) 

(0,0,40) 

(0,0,65) 

(0,0,13.2) 

20* 

0 

l.C.1.3 

10 

100 

1000 

100 

(0,0,10) 

(0,0,40) 

(0,0,65) 

(0,0,13.2) 

40* 

0 

l.C.2.2 

10 

100 

1000 

100 

(0,0,10) 

(0,0,40) 

(0,0,65) 

(0,0,13.2) 

10 

0.5 

l.C.2.4 

10 

100 

1000 

100 

(0,0,10) 

(0,0,40) 

(0,0,65) 

(0,0,13.2) 

10 

1 

l.C.2.5 

10 

100 

1000 

100 

(0,0,10) 

(0,0,40) 

(0,0,65) 

(0,0,13.2) 

10 

1.5 




2. 

In-plane Fields {B x , 

, 0, 0) and (0, B y , 0) 




2.B. Isolated Cloud with Clump 








2.B.1.0 

10 

100 

1000 

- 

(40,0,0) 

(40,0,0) 

(40,0,0) 

- 

- 

- 

2.B.1.1 

10 

100 

1000 

- 

(10,0,0) 

(10,0,0) 

(10,0,0) 

- 

- 

- 

2.B.1.2 

10 

100 

1000 

- 

(65,0,0) 

(65,0,0) 

(65,0,0) 

- 

- 

- 

2.B.2.0 

10 

100 

1000 

- 

(0,40,0) 

(0,40,0) 

(0,40,0) 

- 

- 

- 

2.B.2.1 

10 

100 

1000 

- 

(0,10,0) 

(0,10,0) 

(0,10,0) 

- 

- 

- 

2.B.2.2 

10 

100 

1000 

- 

(0,65,0) 

(0,65,0) 

(0,65,0) 

- 

- 

- 

2.C. Cloud 

Collision 










2.C.1.0 

10 

100 

1000 

100 

(40,0,0) 

(40,0,0) 

(40,0,0) 

(40,0,0) 

10 

0 

2.C.1.1 

10 

100 

1000 

100 

(10,0,0) 

(10,0,0) 

(10,0,0) 

(10,0,0) 

10 

0 

2.C.1.2 

10 

100 

1000 

100 

(65,0,0) 

(65,0,0) 

(65,0,0) 

(65,0,0) 

10 

0 

2.C.2.0 

10 

100 

1000 

100 

(0,40,0) 

(0,40,0) 

(0,40,0) 

(0,40,0) 

10 

0 

2.C.2.1 

10 

100 

1000 

100 

(0,10,0) 

(0,10,0) 

(0,10,0) 

(0,10,0) 

10 

0 

2.C.2.2 

10 

100 

1000 

100 

(0,35,0) 

(0,65,0) 

(0,65,0) 

(0,65,0) 

10 

0 

3. Mixed Fields (B x , B y: B z ) 

3.B. Isolated Cloud with Clump 








3.B.1 

10 

100 

1000 

- 

(10,0,0) 

(10,0,38.7) 

(10,0,64.2) 

- 

- 

- 

3.B.2 

10 

100 

1000 

- 

(0,10,0) 

(0,10,38.7) 

(0,10,64.2) 

- 

- 

- 

3.C. Cloud Collision 










3.C.1.0 

10 

100 

1000 

100 

(10,0,0) 

(10,0,38.7) 

(10,0,64.2) 

(10,0,12.9) 

10 

0 

3.C.1.1 

10 

100 

1000 

100 

(10,0,0) 

(10,0,38.7) 

(10,0,64.2) 

(10,0,12.9) 

5 

0 

3.C.1.2 

10 

100 

1000 

100 

(10,0,0) 

(10,0,38.7) 

(10,0,64.2) 

(10,0,12.9) 

20* 

0 

3.C.1.3 

10 

100 

1000 

100 

(10,0,0) 

(10,0,38.7) 

(10,0,64.2) 

(10,0,12.9) 

40* 

0 

3.C.2.0 

10 

100 

1000 

100 

(0,10,0) 

(0,10,38.7) 

(0,10,64.2) 

(0,10,12.9) 

10 

0 

3.C.2.1 

10 

100 

1000 

100 

(0,10,0) 

(0,10,38.7) 

(0,10,64.2) 

(0,10,12.9) 

5 

0 

3.C.2.2 

10 

100 

1000 

100 

(0,10,0) 

(0,10,38.7) 

(0,10,64.2) 

(0,10,12.9) 

20* 

0 

3.C.2.3 

10 

100 

1000 

100 

(0,10,0) 

(0,10,38.7) 

(0,10,64.2) 

(0,10,12.9) 

40* 

0 

3.D. Cloud Collision with Impact Parameter 







3.D.0 

10 

100 

1000 

100 

(10,0,0) 

(10,0,38.7) 

(10,0,64.2) 

(10,0,12.9) 

10 

0.5 
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and inside the clump. We set the scalar S to 1 inside the 
clump and 0 outside. We added a conservation equation 
in Enzo 2.0 to advect the scalar, i.e., 

+ V.(/>5v) = 0, (1) 

with p the cell density and v the velocity. 

We model a numerical domain of 256 2 pc 2 which is 
covered by a uniform grid of 1024 2 , giving a grid cell 
size of 0.25 pc. For the fiducial model, two additional 
AMR grid levels are included, thus increasing the effec¬ 
tive resolution to 4096 2 with a grid cell size of 0.0625 pc 
on the finest level. This resolution is sufficient to study 
the transition from subcritical to supercritical clouds and 
clumps. 

We use several criteria to determine refinement: a cell 
is refined when there is a strong local gradient of variables 
(i.e., when the relative slope | q(i + 1 ) — q(i — l))/q(i)\ 
across variable q at index i exceeds 0.4), when it is part 
of a shock front (defined by a relative pressure jump of > 
0.33), and when the local Jeans length is not covered by 
at least 4 cells (needed to avoid artificial fragmentation 
[Trnelove et al.|[ l997). 

3. HEATING AND COOLING FUNCTIONS 

We model the thermal properties of the ISM using 
a Photo-Dissociation Region (PDR)-based method that 
follows and expands upon VLBT2013 and is detailed in 
the following subsections. 

3.1. Implementation of Thermal Processes 

Implementation into the Enzo code involves calculating 
the net heating rate for a given cell: 

H = nn[r — tirA] erg cm -3 s -1 , (2) 

where T is the heating rate and A is the cooling rate. The 
net cooling rate introduces a cooling timescale: t coo \ = 
A'int /1 H |. The internal energy of the gas is defined by 
Emt =p /(7 - !)• 

We adopt a mean particle mass of p = 2.33 ran (valid 
for molecular gas with 1 He per 10 H and ignoring con¬ 
tributions from other species). For simplicity, this value 
of p is adopted through the simulation domain, i.e., even 
in the ambient, “atomic” medium. 

The dynamics of the simulation (and the objects of 
main interest) are dominated by gas at densities of tzr > 
10 2 cm -3 , which correspond to equilibrium temperatures 
of - 10 K (details described below). Thus, we adopt the 
value of 7 = 5/3 for the entire simulation domain. While 
this does not account for the excitation of rotational and 
vibrational modes of H 2 that would occur in shocks, we 
consider that this is the most appropriate single-valued 
choice of 7 for our simulation set-up, given our focus on 
the dynamics of the dense molecular gas. 

The chosen values of 7 and p set sound speeds of 
C s = ( 7 fcr//imp ) 1 / 2 « 0.24(T/10 K ) 1 / 2 km s" 1 . Since 
Cool is often shorter than the hydro dynamical time, the 
temperature and internal energy are sub-cycled and up¬ 
dated, assuming constant density, until the hydrodynam- 
ical time step is reached. This is more computationally 
efficient as a method for preventing excessive heating or 
cooling, than evolving all variables on timesteps equal to 
the cooling or heating times. 


3.2. Density-Extinction Relation 

Simulating a ~ kpc 3 * * region of a galactic disk, 
VLBT2013 found a monotonically increasing relation be¬ 
tween density and average (six orthogonal ray) column 
extinctionj_which defined an effective visual extinction 
(Glover & Mac Low 2007). This relation was resolution- 
limited at high densities due to the effective visual ex¬ 
tinction being dominated by absorption within a single 
0.5 pc cell. For the simulations performed in this pa¬ 
per, we use a modified extinction curve normalized to 
estimated values of the Warm Neutr al Med ium (WNM): 
Ay ~ 0.01 mag for nn = 0.03 cm -3 ( Wolfire et al.|2QQ3 ), 
GMCs: Ay ~ 1 mag for tir = 100 cm -;J 7 and starless 
cores: Ay ~ 30 mag for uh = 10 6 cm -3 . 

To fit these constraints as well as retain the physical 
relationships represented in the VLBT2013 curve, we ig¬ 
nore the effects of individual cell extinction at high den¬ 
sities and instead perform a logarithmic extrapolation 
from nn = 10 3 cm -3 . This intermediate curve is fitted 
to the normalization points via a smooth scaling func¬ 
tion, producing the final density-extinction relation (see 

Fig .[§. 



Fig. 2.— Average visual extinction as a function of density. The 
red solid line represents the adopted Ay versus tir relation, based 
on three observational constraints (see text). For comparison, the 
blue dashed line represents the relation used by VLBT2013, with 
the dotted line showing the resolution limit due to extinction within 
the cell itself. 


3.3. Non equilibrium PDR Heating and Cooling Rates 

We next utilize both PyPDR 0 (described below) a nd 
Cloudy (version 13.02, last described in Ferland et al. 
(20l3|)), photoionization simulation codes, to generate 
tables of non-equilibrium heating and cooling rates as 
functions of density, temperature and radiation field in¬ 
tensity. Our default value of FUV radiation field in¬ 
tensity is Go = 4, following conditions developed for 
the ^ 4 kpc molecular ring region of the inner Milky 
Way (VLBT13). Then, given the Ay vs. ur relation 
described in the last subsection, each value of density 
has a unique value of received FUV intensity, allowing a 
2D (ur,T) grid of heating and cooling rates to be suffi- 

1 http://www.mpe.mpg.de/~simonbr/research_pypdr/index.html 
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cient. However, to calculate this 2D grid self-consistently 
does require calculation of arbitrary PDR models with 
input density, temperature and radiation field, in or¬ 
der to calculate species abundances correctly, especially 
of molecules like H 2 and CO that have abundances set 
by FUV photons whose propagation is affected by self¬ 
shielding. 

To carry out these PDR models we primarily utilize 
PyPDR, which is a minimal Python-based PDR code 
that self-consistently calculates chemical, thermal, and 
molecular properties within a slab of gas irradiated with 
FUV photons. 

PyPDR eimplements the same chemical network as in 
|Rollig et al. (2007|), which includes reactions for H 2 for¬ 
mation, cosmic ray induced reactions, photodissociation 
(including self/mutual-shielding of H 2 , CO and C) and 
gas-phase reactions. 

The following heating and coolin g rates are imple¬ 
mented: H 2 pump ing an d li ne coo ling (Rohig et al.|2006), 


H 2 for mation (Sternberg & Dalgamo 


ation (Jonkheid et al 


(Tielens 2005), photoe 


.1989}, H 2 dissoci- 
. 2004}, gas-grain heating/cooling 
electric 


heating and recombi natio n 


cooling ([H akes & Tielens||1994 ), Ly-<a cooling (Sternberg 


fc Dalgarno|l989 ) , optic al Oxygen-6300 A cooling fl Stern- 
berg &; Dalgarno 198 9]), heating by C-ioniz ation (| Black 
V van Di shoeck 1987; Jonkhei d et al.||2004| ), cosmic ray 


heating ([jonkheid et al. 2004), line cooling by 01, CII, 
Cl (fine structure), CO and i3 CO (rotational) calculated 
from the non-LTE excitation of OI, CII, Cl, CO, and 
13 CO using an escape probability approach . Data from 
the LAMDA database ( Schoeier et al.| |2005) is used. 

While the PyPDR chemical network includes only ^30 
atoms and molecules, it still performs well in bench- 
mar ktests^producing results similar to larger PDR codes 
(Rol lig et al.|20 07). However, as PyPDR was developed 
for temperatures only up to ^ 10 4 K, we do not use it 
for higher temperature conditions. 

Cloudy, on the other hand, follows a much larger 
number of species than PyPDR and can treat T > 10 4 K 
gas. Thus, we utilize it in this regime. However, for our 
purposes of defining non-equilibrium heating and cool¬ 
ing functions that utilize a two step process where high 
spectral resolution line self-shielding output is needed as 
a general input for the next PDR calculation, the public 
version of Cloudy does not automatically provide such 
output information. Thus we have adapted the PyPDR 
code of Bruderer for this purpose. 

We set up the density-temperature parameter space for 
the arrays as follows. The density range is nn = 10 -3 
to 10 6 cm -3 , in steps of 0.1 dex (91 values) while the 
temperature spans from 2.7 K to 10 5 K in steps of 0.046 
dex (100 values). PyPDR was used to calculate the bulk 
of the rates, from T = 2.7 to 10 4 K, while Cloudy was 
used for T = 10 4 to 10 5 K. 

The procedure for both PDR codes generally follows 
that of VLBT201T w hich used Cloudy version 8.02 and 
was based off of Smith et al. (2008). Any code-specific 
differences will be mentioned in the relevant sections. 

First, the unextinguished local interstellar radiation 
field (ISRF) with Go = 4 is incident on an absorb¬ 
ing slab of gas with abundances, metallicities, and dust 
resembling that of the local ISM. We include the cos¬ 
mic microwave background radiation as well as a back¬ 


ground of cosmic rays with primary ionization rate of ( = 
1.0 x 10 -16 s -1 . The column density of the gas slab is de¬ 
termined by the previously described density-extinction 
relation and the linear relation between column density 
and visual extinction ( Ay = 5.35 x 10 -22 IVh mag). For a 
given density, a PDR model is calculated through a slab 
with column corresponding to the particular density. In 
our case, the density of the slab also follows this given 
density, as the extinction should be dominated by the 
local density for GMC regions. Note, in VLBT13, the 
density of the slab was fixed at nn = 1 cm -3 . 

At the depth of the specified column (i.e., the final 
cell of the PDR model) the temperature of a parcel of 
gas is then varied, with heating and cooling rates calcu¬ 
lated for the specific temperature. The calculations are 
repeated for the entire temperature range, yielding the 
temperature and density dependent heating and cooling 
functions. 

In PyPDR, the code allows self-consistent calcula¬ 
tion of general, nonequilibrium heating and cooling rates 
given species abundances set by equilbrium PDR condi¬ 
tions for given extinction, density and equilibrium tem¬ 
perature. The key PyPDR results for T eq , H 2 , and CO 
abundances are shown in Fig. [3] 

Arrays of heating and cooling rates were created using 
both PyPDR (T = 2.7 K to 10 4 K) and Cloudy (10 4 K 
to 10 5 K), then smoothly joined along the temperature 
dimension using the function: 

R(T) = 10.0 loSlo ^ c ^ T ^ a ^ T ^ +1 ° Slo ^ jRp ^ T ^ 1,0_a ^ T ' ) ^ (3) 

where R(T) is the final, smoothly combined rate, calcu¬ 
lated from T the gas temperature, Rc the Cloudy rate, 
Rp the PyPDR rate, and the Fermi-Dirac smoothing 
function: 


a(T) = 


1 + exp [—10.0(log 10 T - 4.0)] 


(4) 


This joined the functions at T = 10 4 K, where good 
agreement still occurred between the models. 

From the resulting final arrays of heating rates and 
cooling rates, a bilinear interpolation is performed to de¬ 
rive rates for any density and temperature combination. 
The final, combined 2D interpolation plots for cooling, 
heating, and net heating as functions of density and tem¬ 
perature are displayed in Fig. [4] These plots show the 
total cooling, heating, and net rate at all densities and 
temperatures. 

Note, that from T = 10 5 K up to T = 10 8 K, we ignore 
heating and a dopt the cooling rates derived by |Sarazin fc| 
White (1987). Values beyond the array domain adopt the 
limiting values. Thus, for any cell in our Enzo simulation, 
the density and temperature are read in and cooling and 
heating rates are returned. 


3.4. Heating and Cooling Components 

A breakdown of specific heating and cooling compo¬ 
nents at the equilibrium temperature is shown in Fig. [5] 
Photoelectric heating of dust grains is the dominant heat¬ 
ing source for low-density gas up until nn ^ 10 2 cm -3 . 
Above this density, the higher dust extinction blocks ex¬ 
ternal FUV photons and thus reduces photoelectric heat¬ 
ing. The ubiquitous flux of cosmic rays then becomes the 
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Fig. 3.— (top) PyPDR equilibrium temperature as a function of 
density. The density cor resp onds to a value of Ay from from Fig. [2] 
Details are discussed in §3.3| Lines of constant pressure are plotted 
in gray to more easily show regions of thermal instability, (middle) 
H 2 fraction as a function of density. Hydrogen becomes essentially 
fully molecular at densities above uh — 80 cm -3 , (bottom) CO 
fraction as a function of density. The carbon becomes fully molec¬ 
ular in the form of CO at densities above wh ~2x 10 3 cm -3 . 

main heating component in high-density gas. H 2 forma¬ 
tion also contributes as the cloud becomes fully molecu¬ 
lar. 

The main coolants in the low density, ionized/atomic 
region (njj < 1cm -3 ) are Ly-<a and Hydrogen recombina¬ 
tion lines. As density increases, various atomic lines (01, 
CII, Cl) become dominant coolants. These species inelas- 
tically collide with H and He, exciting internal degrees 
of freedom and subsequently decaying through photon 
emission. Molecular lines (CO, 13 CO) provide large con¬ 
tributions in cooling as density increases, temperature 
decreases, and the gas reaches high levels of molecular 
abundance. At the highest densities (> 10 4 cm -3 ), gas- 
grain cooling dominates as collisions between dust grains 
and gas molecules lead to emission of infrared photons 
from the decay of lattice vibrations. 



n„ [cm 3 ] 


Fig. 4.— Density and temperature dependent interpolated ar¬ 
rays of (top) cooling, (middle) heating, and (bottom) ratio of heat¬ 
ing/cooling. The contours show constant rates (top and middle 
panels) and ratios (bottom panel), e.g., in the ratio map, the 10° 
contour represents the equilibrium temperature. 

3.5. Observational Diagnostics 

In addition to providing a better understanding of the 
dominant physical processes occurring at different den¬ 
sities and temperatures, the heating/cooling component 
breakdown also enables the creation of observational di- 
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Fig. 5. — Component breakdown of the main cooling (top) and 
heating ( bottom ) rates per unit volume as a function of density at 
the respective equilibrium temperatures (given in Fig. [ 3 J. 


agnostics in the form of line emissivities. Here we focus 
on a preliminary investigation into high-J CO to see if 
they are good diagnostics of shocks arising from cloud 
collisions. The PDR-derived cooling data include contri¬ 
butions from the first 40 rotational lines for both 12 CO 
and 13 CO. Similar to the method of creating the cooling 
and heating functions, tables of density and temperature 
dependent emissivities were compiled to allow calcula¬ 
tion of observable quantities in the form of integrated 
intensity maps and spectra. 

Via post-processing, integrated intensity maps can be 
derived from the volume emissivity functions coupled 
with simulation outputs. We assume a fixed 1 pc thick¬ 
ness of the simulation volume for calcuation of these 
maps. Given this simplistic, highly-idealized 2D geom¬ 
etry of cloud structures presented in this initial paper, 
for simplicity we do not calculate radiative transfer of 
the emissivities from each cell, but simply sum their 
contributions as if their emission reached us with neg¬ 
ligible attenuation. However, note that the emissivities 
of lines from PyPDR do already account for an ideal¬ 
ized cloud optical depth via an esca pe p robability for¬ 
malism through the PDR layer (see § 3.3). Further de¬ 
tailed study of observational diagnostics of cloud-cloud 
collisions based on 3D simulations and including full ra¬ 
diative transfer will be deferred to a future paper. 


For ease of comparison with Galactic clouds, we as¬ 
sume a fiducial cloud distance of d = 3 kpc and depth 
of 2 = 1 pc. We use the line volume emissivities derived 
from the PyPDR to determine an integrated intensity 
for each cell in the simulation. Integrated intensities are 
derived via: 

2k C 

= ^ J T mb diA (5) 

where I v is the specific intensity, A is the wavelength of 
the chosen molecular line, and T mb is the main beam 
temperature. Changing variables from v to v and sub¬ 
stituting, we have 

/ r - d ” = S / =8Sr (6) 



where j is the volume emissivity, V is the cell volume, 
and Cl is the solid angle subtended by the cell. 

While values of z and d are assumed in order to provide 
some observational outputs, the intensity maps can be 
scaled for any desired thickness, given the optically-thin 
assumption. Integrated intensity maps of CO lines and 
line ratios with rotational excitations J=2-l, 3-2, and 8-7 
using this method are presented and discussed in m 

In addition to integrated intensity maps, spectra of 
the corresponding observational volumes can be created, 
simply by plotting the distribution of specific intensity as 
a function of line of sight velocity. Synthetic spectra of 
the 13 CO(J=2-l), 13 CO(J=3-2), and 12 CO(J=8-7) lines 
for an isolated GMC and a GMC collision case, viewed 
along sight lines within the 2D simulation plane, are com¬ 
pared. Velocity gradients derived from these spectra are 
described as well in Sj5] 

4. RESULTS 


4.1. Out-of-plane magnetic fields 

Here we assume the magnetic fields are orientated or 
thogonal to the 2D simulation plane and t hus the coll i 
sion veloc i ty. U sing the virial theorem, Chandrasekhar 


& Fermi (1953) showed that magnetic fields support 


the cloud preventing gravitational collapse if the aver¬ 
age magnetic field strength exceeds 


B crit = 2nRpG 1/2 , 


(7) 


where p is the average density of the cloud. Note this 
assumes the external magnetic field to be negligible. We 
systematically vary the magnetic field strength to probe 
the sub- and supercritical regimes and to understand the 
transition from sub- to supercritical. 


4.1.1. Isolated cloud 

The simplest case to consider is an isolated cloud. 
For our adopted parameters (see GMC 1 values of Ta¬ 
ble IT] and runs l.A.x in Tab. [2|, the critical magnetic 
fielais 27.9 pG. We initialize the cloud with a uniform 
out-of-plane magnetic field, sampling values of 10, 20, 
27.9 and 40 pG and setting the ambient field to 10 pG. 
We also carry out an unmagnetized simulation. Cloud 
evolution is followed for 10 Myr, more than 2 freefall 
times, tff = (37r/[32Gp]) 1 / 2 , which is ~ 4.35 Myr for the 
adopted cloud values. Note that this expression for the 
free-fall time is for a uniform sphere and not for an in¬ 
finite, uniform cylinder. This is an intentional choice to 
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Fig. 6 .— Average density of the cloud over time for initial B 
strengths of 0, 10, 20, 28, and 40 /iG. The blue dash-dotted line 
marks when the B = 0 /jlG case is affected by numerical effects. 
Critical field strength occurs at B C rit = 28/iG. The dotted vertical 
line denotes the freefall time of the cloud. The dotted horizontal 
lines show p cr it predicted for each model. 


use a common definition, as we will extend this work to 
3D in future studies. 

Figure [6] shows the evolution of the average cloud den¬ 
sity for each mod el, tracked using the advective scalar 
method (see j |2.2| ). 

For the unmagnetized, pure hydro dynamical model, 
the line-mass of the cloud exceeds the critical value 
2 a 2 /G and thus collapses unimpeded. All of the cloud 
mass ends up into a single cell in about a free-fall time 
and then continues to slowly accrete mass from the ex¬ 
ternal medium (as seen in the shallow flat slope of the 
mean density at 5 Myr). Note that after t ~ 6Myr, the 
cloud material is no longer tracked well, likely the result 
of a numerical artefact due to numerical diffusion. Note 
that we track the cloud (or clump, as in later cases) by 
advecting a scalar field S which we set to 1 inside the 
cloud and 0 outside. The cloud is defined by material 
with S > 0.5. By the time the cloud collapses, most of 
the mass is within a small number of grid cells. As evolu¬ 
tion continues, cloud and ambient material mixes so that 
the scalar is now below the defined value and its mass is 
no longer accounted for. 

If magnetic fields are present, the magnetic pressure 
supports the cloud against gravitational collapse. The 
oscillation is due to the transition towards an equilibrium 
density and magnetic field distribution. 

For magnetic field values below the critical value, the 
cloud is initially supercritical and thus collapses. The av¬ 
erage density follows the same evolution as for the pure 
hydrodynamic model. However, the clouds do not col¬ 
lapse completely, even for field strengths close to the ex¬ 
ternal magnetic field value. Actually, as long as a mag¬ 
netic field is present, the collapse of the cloud is impeded. 
This can be easily understood from conservation of mag¬ 
netic flux and mass in a 2D geometry. Both the average 
magnetic field and density are proportional to 1/R 2 . The 
critical magnetic field, however, is proportional to 1/R 
as F> cri t oc Rp (see Eq. \7l. Thus, as the cloud collapses, 
the average magnetic field in the cloud increases faster 


than the critical magnetic value. While the initial mag¬ 
netic field is too weak to support the cloud, subsequent 
contraction causes the field to eventually become strong 
enough to halt collapse. Thus, the cloud transitions from 
a supercritical regime to a subcritical one. 

Using the initial line mass and magnetic flux of the 
cloud, the mean density at which the transition from 
supercritical to subcritical occurs, is given by 


Pcrit — 


47t mfG 

<p ' 


( 8 ) 


For B = 10 pG and nn = 100 cm -3 , we find n cr it = 
776 cm -3 . Figure M indeed shows that the gravitational 
collapse oscillates close to this value for B = lOpG, grad¬ 
ually settling towards an equilibrium state. The final 
densities are slightly higher than predicted, presumably 
due to external pressure from the ambient, magnetized, 
infalling gas. 

This result is specific to the 2D cylindrical geometry 
adopted here. For a sphe rical cl oud, th e critical ma gnetic 
field strength is given by [Chandrasekhar & Fermi ( 1953| ) 


B c 


2.5tt RpG 1/2 . 


(9) 


Note the similarity with the expression of the critical 
value for a cylindrical cloud (eq. [7]) . If the cloud is ini¬ 
tially supercritical, we can assume nearly isotropic col¬ 
lapse. Then, the density is proportional to 1/R 3 , so that 
H cr i t oc 1/R 2 . Because of magnetic flux freezing in ideal 
MHD, the mean magnetic field of the cloud is proporp- 
tional to 1/R 2 . This means that the critical magnetic 
field and the mean magnetic field have the same propor- 
tianlity and the cloud remains supercritical during the 
collapse. 

It is clear that, if high gravitational collapse is to be 
achieved in 2D, flow along field lines (along the direction 
of collapse) or flow through field lines (due to, e.g., am- 
bipolar diffusion or turbulent reconnection) must occur. 
Alternate field geometries are discussed in later sections. 


4.1.2. Isolated cloud with embedded clump 

We now embed a clump within the cloud discussed in 
the previous section. The critical magnetic field for the 
clump is 65 pG, while the critical value for the cloud 
has increased to 40 pG (as the average density of the 
cloud is higher with the embedded clump). We exam¬ 
ine the effect of various magnetic field strengths in the 
cloud and clump. First, we keep the the clump mag¬ 
netic field constant and vary the cloud value. This tells 
us more about the evolution of an equilibrium clump in 
a sub- or supercritical cloud. Then, we keep the cloud 
magnetic field constant while varying the clump value. 
These correspond to runs l.B.x in Tab. [2j Although we 
are restricted with our cloud and magnetic field geome¬ 
try, these results are useful for understanding more com¬ 
plex simulations. Results are shown in Figures [7| and 

E 

For a constant Bqmc near its critical value, the evolu¬ 
tion of the clump is entirely determined by the ratio of 
its gravitational and magnetic energy (its thermal energy 
is negligible). Using eq.j8l we find that, for B c \ = 40pG, 
the average density of the clump increases by a factor of 
2.7. However, Figure [7] shows an increase twice this value. 
This is due to the initial contraction of the cloud between 
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Fig. 7. — Average clump density (top panel) and temperature 
(bottom panel) versus time for constant GMC magnetic field B\ = 
40 /aG and varying B c \ = 40, 55, 65, 75 ^G. Critical field strength 
occurs at B c \ = 65 gG. The solid line represents the average 
value in the clump defined by the scalar value S > 0.5, while the 
shaded regions show the averages between S > 0.25 and S > 0.75. 
This convention for clump definition is followed throughout the 
remainder of the paper. 



Fig. 8. — Average clump density (top panel) and temperature 
(bottom panel) over time for constant clump magnetic field B c \ = 
65 gG and varying GMC field strength as B\ = 10, 20, 40, 65 gG. 
Critical field strength occurs at B\ = 40 gG. 


0-2 Myr as it tries to set up an equilibrium distribution. 
After this initial adjustment phase, the average density 
of the clump drops to a few times the initial value as 
expected. For higher initial magnetic fields in the clump 
the density increases by a smaller factor, as also expected 
from eq. [8| 


For a constant B c i, changes in the average clump den¬ 
sity are driven by external pressure from the surrounding 
GMC material. This external pressure is relatively larger 
for the lower GMC magnetic fields. In these cases, the 
GMC is initially supercritical and starts to contract grav¬ 
itationally. The average density of the clump increases 
maximally by a factor of ^ 5. The clump magnetic field 
is strong enough to resist gravity but the clump is fur¬ 
ther compressed to higher densities because the external 
pressure contributes non-negligibly to the gravity. For 
stronger GMC magnetic fields, however, the density of 
the clump is not increasing because of the pressure ex¬ 
erted by the GMC. Instead, the clump is initially no 
longer subcritical. The external (i.e., GMC) magnetic 
field is not negligible and should be taken into account 
when deriving the critical value. For high GMC magnetic 
fields, the critical value of the clump is actually greater 
than 65 /iG, and thus it initially collapses gravitationally. 
However, it is not highly supercritical so the density in¬ 
crease is quite modest. At the same time, GMCs with 
higher magnetic fields expand after 2-2.5 Myr (see previ¬ 
ous section). The external magnetic field then decreases, 
as well as the critical magnetic field of the clump. This 
results in re-expansion of the clump to near its initial 
value. Our results suggest that increasing external pres¬ 
sures is a possible way to trigger a sub-to-supercritical 
transition. 


4.1.3. Colliding clouds: head-on collisions 

A significant source of additional pressure can be pro¬ 
vided by ram pressure of cloud collisions and the result¬ 
ing thermal and magnetic pressure released in shocks. 
The ram pressure depends on the relative collision speed, 
u re i 2 . We investigate different collision speeds, i.e., rv e i = 
5,10, 20, and 40 km s -1 (see runs l.C.l.x in Table [ 2 ]) . 

Density and temperature slices at different stages of 
the evolution are shown for v re \ = lOkms -1 in Figure [9j 
The two clouds are initially separated such that the col¬ 
lision occurs at 4 Myr, which allows for an initial re¬ 
distribution of the density in the cloud (see Figures 10 
and 12). The cloud-cloud collision compresses the clouds 
and the clump, leading to higher densities. The collision 
also gives rise to many shocks propagating through the 
clouds. Such shocks contribute to raising the pressure 
around the clump. High-temperature shock fronts are 
present within the otherwise cold (~ 15 K) clouds. The 
magnitude of the magnetic field also increases as material 
is compressed. This increase in magnetic pressure pre¬ 
vents the clump from collapsing completely, even with 
the additional external pressure of the collision. 

Figure [lO] shows the evolution of the average clump 
density and temperature for different collision speeds 
compared to the non-colliding case. Note that, for 
u re i = 20 and 40 km s -1 , the collisions were performed 
in the reference frame of the clump, to avoid high flow 
velocity induced numerical diffusion effects that can have 
modest effects on clump boundary definitions, mostly af¬ 
fecting measurement of clump temperature. 

Due to the utilized set-up, a collision front between 
the low-density ambient envelopes arises in between the 
clouds. Before the clouds directly interact, they are be¬ 
ing influenced somewhat by this high pressure post-shock 
collision region. However, the pressure here is much less 
than the ram pressure resulting from the GMC-GMC col- 
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Fig. 9. — Evolution of clouds colliding head-on (zero impact parameter) with snapshots shown at 2.05, 4.01, 4.99, 5.96, and 7.92 Myr (see 
run 1.C.0 in Table [2]) Here, magnetic fields are near critical values and directed out-of-plane. B c \ = 65 g, G, B\ = 40 fiG , and Bo = 10 g G. 
(top row ) Maps of wh and (bottom row) temperature, with black vectors representing velocity are shown. The advective scalar defining the 
clump is shown by grey contour lines, representing the scalar value S = 0.25, 0.5, and 0.75. 


lision, given the factor of 10 difference in GMC to ambi¬ 
ent density. The effects of the shocked ambient medium 
on the clump can be seen in Fig. 10, at t ;$ 4 Myr. 
Both density and temperature are affected, more notice¬ 
ably at 20 and 40 km s _1 , but the ensuing GMC collision 
dominates the subsequent clump evolution. These effects 
due to the shocked ambient medium can be seen in the 
density and temperature evolution for colliding cases in 
subsequent runs, discussed below. 

As expected, the density and temperature of the clump 
resulting from the GMC collision increase with collision 
speed, as higher velocities induce stronger shocks with 
larger compressions. However, even for v re \ = 40 km s -1 , 
the increases in clump density are only modest: about 
a factor of 2 to 3 times greater than the isolated case. 
Of course, some of this is due to the specific geometry 
we adopt here. For other cloud geometries, e.g., spher¬ 
ical clouds in 3D, and magnetic field geometries, e.g., 
more parallel to collision velocities, this extra pressure 
may yet be sufficient to trigger the transition from sub- 
to supercritical. The collision models do show larger ex¬ 
cursions in clump mean temperatures, which would be 
expected to have an impact on astrochemical processes 
in the clump. 


4.1.4. Colliding clouds: off-axis collisions 

Off-axis collisions, in which the impact parameter was 
varied, were also explored. GMC 2 was placed at dif¬ 
ferent perpendicular distances, 6, from GMC l’s line of 
symmetry and the mean density of the clump material 
was tracked over 10 Myr. Figure |TT] displays the mor¬ 
phology of the collision for b = 0.5i?i. The clouds inter¬ 
act at ^ 4 Myr in an asymmetric manner, creating fil¬ 
amentary structures and imparting angular momentum 
on the coalesced structure. Compared with the on-axis 
head-on collisions, the resulting structures are morpho¬ 
logically more filamentary but the level of gravitational 
contraction is roughly equivalent. The lack of complete 
gravitational contraction is expected because of the flux¬ 


freezing limitation of out-of-plane fields described above. 
In addition, any angular momentum in the final struc¬ 
ture also helps to support the clump, further reducing 
the degree of contraction. 

The average clump densities for various impact param¬ 
eters are compared in Figure [T2| Collisions at ^ 4 Myr 
show varying factors of density increase, with higher av¬ 
erage densities for smaller values of b (more direct col¬ 
lisions). As with the case of head on collisions, clump 
densities are only increased by at most a factor of a few. 


4.2. In-plane magnetic fields 

The primary inhibitor of complete collapse in the out- 
of-plane (B z ) magnetic field runs is flux freezing, i.e., 
gas cannot move across magnetic field lines. Therefore, 
in this section, we change the direction of the magnetic 
field from orthogonal to the plane to be within the plane. 
Contrary to the out-of-plane field models where the mag¬ 
netic field value is higher inside the cloud than outside it, 
we assume a uniform magnetic field across the cloud and 
external medium. For such clouds, gravitational collapse 
proceeds preferentially along the magnetic field lines. 
While forces supporting the cloud are much great er per¬ 
pendicular to the magnetic field lines, Tomisaka (2014) 
showed that a uniform in-plane field geometry can yield 
magnetically subcritical configurations for infinite cylin¬ 
ders. The maximum line mass was evaluated as 


Vax — 22.4 


Re) 


0.5 pc ) V10/iG 


Bo 


Mq pc 


-1 


+13.9 


( -^ M® pC 1 (10) 

VO.19 km s -1 / OF v J 


for the isothermal case. 

Multiple field strengths were explored, sampling val¬ 
ues previously used in the out-of-plane cases to keep the 
total magnetic pressure component consistent. We ap¬ 
ply \B\ = 10,40, and 65 fiG and analyze the effects on 
isolated and colliding cases. These models correspond to 
runs 2.x in Table [2 













12 


Wu et al. 



0 2 4 6 8 10 


Time [Myr] 

Fig. 10.— Average clump density (Top panel) and temperature 
(2nd panel ) over time comparing the effects of collision velocity for 
out-of-plane field geometries (see runs l.C.l.x). Here, out-of-plane 
magnetic fields are near critical (B c \ = 65 gG, B\ = 40 g G, and 
Bo = 10 gG). Collision velocities u re i = 5,10, 20, 40km s _1 are 
shown, along with the evolution of the isolated GMC with clump. 
The ratios of the colliding cases compared to the isolated case for 
average density (3rd panel) and temperature (bottom panel) are 
also shown. 


4.2.1. Isolated cloud with embedded clump 

With uniform B x and B y fields, the GMC and clump 
collapse along the direction of the field to form dense 
sheets perpendicular to the field lines. The timescales 
associated with their contraction are of the order of the 
spherical free-fall time tff, i.e., ~ 1.6 Myr for the clump 
and 4.4 Myr for the GMC. After the initial collapse par¬ 
allel to the magnetic field, the gas starts to contract per¬ 
pendicularly. Complete collapse of the clump takes much 
longer as the gas motions are perpendicular to the mag¬ 
netic field. 

The isolated case is most similar to the models of 


Tomisaka (2014). However, the embedded overdense 
clump dominates the gravitational collapse of the cloud 


The line mass of the clump is 3450 M 0 pc _1 . Equation 10 
yields A max « 1660 M 0 pc -1 for 65 /iG, and even smaller 
values for 10 and 40 /iG. Thus, the maximum supported 
line mass by in-plane magnetic fields is exceeded, and our 
simulations agree with these results. A field strength of 
~ 135 /iG could be used to support the clump, but this 
case was not explored. 


4.2.2. Colliding clouds 

For colliding clouds we again adopt a fiducial relative 
velocity of 10 km s -1 and study two different in-plane 
magnetic field directions, i.e., parallel (i.e., B x ) and per¬ 
pendicular (i.e., By) to the converging flow. Here we 
set the collision time at t = 0 Myr as there is no more 
stable state to be reached. Further, we only study a sin¬ 
gle collision speed as the dynamics are dominated by the 
gravitational collapse of the clouds and clump. Similar to 
the isolated model, the line mass of the clump and clouds 
exceeds the maximum supportable by thermal and mag¬ 
netic pressures. The clouds collapse into flattened sheets 
perpendicular to the magnetic field lines. The timescales 
of gravitational collapse are again of the order of the free- 
fall time. In these highly collapsed scenarios, the clump 
is no longer well tracked at late times due to numerical 
effects. 


4.3. Mixed field geometries 

While the previous sections describe two extremes, i.e., 
either the cloud is maximally supported by magnetic 
fields (out-of-plane magnetic field) or minimally (in-plane 
magnetic field), we now investigate a combination of the 
two geometries. It represents a more realistic situation as 
expected in 3D, where the magnetic field provides some 
support against gravitational collapse, but cannot halt it 
completely, if the cloud is supercritical. 

In these cases, we assume a uniform in-plane magnetic 
field strength of 10 /iG (along the x-axis (Fig. |~L3| ) or 
the y -axis (Fig. fl4|)). The out-of-plane components are 
chosen such that the total field strength has a magni¬ 
tude equal to its critical field strength (see Table [l] and 
runs 3.x in Table [2]). The external medium has zero out- 
of-plane magnetic field component, preserving the total 
field strength of \Bq\ = 10 /iG. Suc h a field is both 
density-dependent (as observed by Crutcher (2012)) and 
is divergence-free. 

4.3.1. Isolated cloud with embedded clump 

As the out-of-plane magnetic field is near-critical and 
strong enough to stabilize the GMC and clump, the early 
evolution is similar to the out-of-plane case (see Fig¬ 
ure [10). A density (and magnetic) gradient is quickly 


established to form an equilibrium structure. However, 
gas also flows along the in-plane magnetic field. Then the 
line mass of the cloud increases while the magnetic flux 
remains constant. The clump and GMC gradually con¬ 
tract, although the associated timescale is much longer 
than the free-fall time. For a larger ratio of in-plane to 
out-of-plane magnetic field, the evolution is faster as the 
out-of-plane magnetic field is less dominant. 
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Fig. 11.— Time evolution of colliding clouds at 2.05, 4.01, 4.99, 5.96, and 7.92 Myr with b = 0.5R\ (see run l.C.2.2 in Table [5]) Here, 
out-of-plane magnetic fields are near critical (B c \ = 65 p-G, B\ = 40 g G, and Bo = 10 gG). (top row) Maps of nn and (bottom )"maps of 
temperature with black vectors representing velocity are shown. The advective scalar defining the clump is shown by grey contour lines, 
representing the scalar value S = 0.25, 0.5, and 0.75. 


4.3.2. Colliding clouds: head-on collisions 


We perform simulations of clouds colliding in this 
mixed-field geometry, again investigating the effect of 
collision speed (y re \ = 5,10, 20, and 40 km s _1 ). Addi¬ 
tionally, the direction of the in-plane component of the 
magnetic field is varied with respect to collision velocity 
(see runs 3.C.1.X and 3.C.2.X in Tab. [2]) 

Similar to the mixed-field isolated cloud case, the col¬ 
liding clouds threaded by a mix of out-of-plane and in¬ 
plane magnetic fields experience a larger compression 
compared to the purely out-of-plane clouds. However, 
the direction of the uniform component of the magnetic 
field affects the late-time behavior of the clump. 

For partial fields parallel to the collision velocity (i.e., 
x-direction), there are temporary density increases of a 
factor of ^ 2 — 3 during the collision, but the average 
clump density actually decreases slightly at ~ 5 Myr and 
beyond, relative to the isolated case (see Figure 15). This 
rebound effect is greater for higher velocities. T he shocks 
initially compress the clump, then subsequently distort 
it, forming a sickle-like shape. From Figure [l3j we see 
that the original clump is broken apart due to the colli¬ 
sion. The temperature of the clump material is affected 
more significantly as high velocity shocks dominate the 
mean clump temperature, temporarily raising it to ^few 
hundred K. The material cools to ^tens of K in the af¬ 
termath of the collision. 

For partial fields perpendicular to the collision veloc¬ 
ity (i.e., ^/-direction), the behavior is nearly identical for 
pre-collision times t < 4 Myr. However, the different 
magnetic field geometry causes the clump to be com¬ 
pressed in a different manner (see Figures 14 and 16). 
In this case, the collision induces no sickle-shaped struc¬ 
ture, but rather the clump stays relatively compact, with 
the average density increasing, but not rebounding. The 
material in the collisional flow interface region freely falls 
into the overdense remnants of the cloud and clump due 
to the orientation of the B-field. Shocks are continually 
created as the global flow and infalling material interact, 


regulated by the magnetic fields. Late time behavior 
after the collision reveals continuously increasing clump 
densities due to infall, with elevated (T ~ 50 — 100 K) 
but roughly level temperatures. 

4.3.3. Colliding clouds: off-axis collisions 

Our final model is a cloud-cloud collision in the mixed- 
field geometry with an in-plane uniform field of B x = 
10 fiG. We have v re \ = 10 km s -1 and additionally apply 
b = 0.5.Ri to GMC 2 (see 3.D.0 in Tab. |J). 

We designate this as our “fiducial case” and run the 
standard resolution, along with one and two additional 
levels of AMR, giving a maximum effective resolution 
of 0.0625 pc. We compare the effects of different res¬ 
olutions in Figure [l7j Pre-collision densities are quite 
well converged, but begin to deviate as the shock waves 
and clump compression are realized at different resolu¬ 
tions. Larger initial differences are seen in the temper¬ 
atures, where higher resolutions lead to generally lower 
average clump temperatures. This is likely due to the 
initial shock created at the boundaries of the uniform 
clump as the density gradient is established. At higher 
resolutions, the post-shock region contributes less to the 
overall clump material. Additionally, inspection of clump 
contours at the various resolutions revealed slightly dif¬ 
ferent clump boundaries arising from the collision. This 
could partially account for the greater discrepancies at 
later times. While these resolution effects are not in¬ 
significant, the key results - relative changes of a cloud 
collision with respect to the isolated case - retain good 
agreement throughout the majority of the simulation (up 
to ~ 8 Myr). 

Figure [18] summarizes the entire fiducial run, show¬ 
ing time evolution for maps of density, tempera¬ 
ture, and common observational bands of 13 CO(J=2- 
1), 13 CO(J=3-2), and 12 CO(J=8-7) as well as the 
12 CO(J=8-7)/ 13 CO(J=2-l) line ratio. These integrated 
intensity maps, based on outputs from the PDR model¬ 
ing as potential observational diagnostics, are discussed 
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Fig. 12. — Average clump density (Top panel) and temperature 
( 2nd panel) over time comparing the effects of impact parameter for 
out-of-plane field geometries (see runs l.C.2.x). Here, out-of-plane 
magnetic fields are near critical (B c \ = 65 pG, B\ = 40 pG : and 
Bo = 10 pG). Impact parameters of b = 0.5,1.0, and 1.5Ri were 
explored. The ratios of the colliding cases compared to the isolated 
case for average density ( 3rd panel and temperature (bottom panel) 
are also shown. 


in ^5] 

Broadly speaking, the effect of a finite impact param¬ 
eter for the GMC collision results in a shearing velocity 
field and asymmetric morphologies as various areas of 
the clump are compressed and distorted. GMC 2 can 
be seen contracting gravitationally as it approaches the 
more massive GMC 1. Prior to the collision, parts of 
GMC 1 and the clump are slightly compressed by the 
bounding shocks arising from the colliding region of the 
ambient material. The collision itself compresses parts of 
the clouds even further, as GMC 2 enters GMC 1 and im¬ 
pacts the clump from the north. From the density and 
temperature maps, shocks can be seen permeating the 


cloud material and passing through the clump through¬ 
out the entire collision process. At later times, the orig¬ 
inal clump material is distorted greatly and even breaks 
apart into a few pieces, but the densest material remains 
inside the main clump region. 

Peak compression due to the collision occurs near 5 
Myr (third column in Figure 18). We investigate this 
further by zooming in on the clump at this timestep and 
mapping various key quantities. This is shown in Fig¬ 
ure [l9j The density, temperature, magnetic fields, and 
velocity gradient in the regions surrounding and includ¬ 
ing the clump are analyzed. Integrated intensity maps 
are also shown in this figure and discussed in © 

The clump, initially a uniform cylinder, remains rela¬ 
tively distinct and contiguous, though at this timestep 
it is undergoing compression and distortion due to the 
cloud collision. What was once GMC 2 can be seen as 
the denser (few xlO 3 cm -3 ) material that has punched 
into GMC 1 and is impacting the clump from the north. 
The average clump density is n h ^ 10 4 cm -3 , embedded 
in GMC material of ~ 10 2 -10 3 cm -3 . 

The clump temperature, on the other hand, is not par¬ 
ticularly distinct from the surrounding material, gener¬ 
ally at a few 10s of K. Shocks of a few 100s of K are seen 
propagating through the clump and cloud. The high tem¬ 
perature material (^ 10 4 K) due to the strong shock cre¬ 
ated by the collision with GMC 2 has penetrated GMC 
1, but has not reached the clump. 

The magnetic fields can be seen corresponding closely 
to the density morphology of the GMC, with field 
strength generally increasing with density. The B-fields 
have strengths of ~ 100 /rG in the compressed GMC ma¬ 
terial and peak at a few hundred /iG within the clump 
and nearby regions. The initial in-plane fields, uniform 
and directed along the collision axis, remain mostly uni¬ 
form, except for where the GMCs have been disrupted. 
Complex field structures arise within the clump and 
cloud material. In this case, there is a loose correla¬ 
tion between magnetic field direction and the direction 
of infalling gas flow to the clump. 

The velocity gradient map shows detailed structure of 
the many shocks propagating throughout the cloud. The 
strongest gradients can be seen corresponding with the 
shocked GMC-envelope interface, as well as the GMC- 
GMC collision region. The velocity magnitudes show 
some turbulent motion being produced by the collision. 

To illustrate the effects of our treatment of nonequi¬ 
librium heating and cooling, Figure [20] compares the 
differences in temperatures between the nonequilibrium 
cooling/heating functions developed in this paper and 
a cooling/heating curve that assumes equilibrium tem¬ 
peratures. Differences primarily occur in the shocked 
regions, as material is shock heated out of equilibrium. 
The temperature maps, upon which the observational 
diagnostics heavily depend, would exhibit very different 
behavior had only a simple equilibrium cooling/heating 
curve been used. 


5. OBSERVATIONAL DIAGNOSTICS 

Here we briefly outline two potential methods of obser- 
vationally diagnosing GMC collisions, based on emission 
of high- J CO lines. However, given the idealized 2D na¬ 
ture of the simulations presented so far, we defer more 
detailed discussion to a future paper that will consider 
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Fig. 13.— Time evolution of colliding GMCs at 2.05, 4.01, 4.99, 5.96, and 7.92 Myr for the x-directed mixed field geometry (see run 
3.C.1.0 in Table [ 2 ]) Here, the total B-field magnitude is near critical while an additional in-plane uniform field of B x = 10 fiG is applied 
throughout the simulation. (top row) Maps of tt-h with magnetic fields represented by streamlines and (bottom row) maps of temperature 
with black vectors representing velocity are shown. The advective scalar defining the clump is shown by grey contour lines, representing 
the scalar value S = 0.25, 0.5, and 0.75. 
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Fig. 14.— Time evolution of colliding GMCs at 2.05, 4.01, 4.99, 5.96, and 7.92 Myr for the y-directed mixed field geometry (see run 
3.C.2.0 in Table [ 2 ]) Here, the total B-field magnitude is near critical while an additional in-plane uniform field of B y = 10 fiG is applied 
throughout the simulation. (top row) Maps of wh with magnetic fields represented by streamlines and (bottom row) maps of temperature 
with black vectors representing velocity are shown. The advective scalar defining the clump is shown by grey contour lines, representing 
the scalar value S = 0.25, 0.5, and 0.75. 
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Fig. 15.— Average clump density ( Top panel) and temperature 
( 2nd panel) over time comparing the effects of collision velocity for 
the B m i x field geometry case (see runs 3.C.l.x). Here, out-of-plane 
magnetic fields are near critical (. B c \ = 65 pG, B\ — 40 pG, and 
Bo = 10 pG) while an additional in-plane uniform field of B x = 
10 pG is applied throughout the simulation. Collision velocities 
u re i = 5,10,20,40km s _1 are shown, along with the evolution of 
the isolated GMC with clump. The ratios of the colliding cases 
compared to the isolated case for average density ( 3rd panel) and 
temperature ( bottom panel) are also shown. 
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Fig. 16.— Average clump density ( Top panel) and temperature 
( 2nd panel) over time comparing the effects of collision velocity for 
the B m i x field geometry case (see runs 3.C.2.x). Here, out-of-plane 
magnetic fields are near critical (. B c \ = 65 pG, B\ — 40 pG , and 
Bo = 10 pG) while an additional in-plane uniform field of B y = 
10 pG is applied throughout the simulation. Collision velocities 
n re i = 5,10,20,40km s -1 are shown, along with the evolution of 
the isolated GMC with clump. The ratios of the colliding cases 
compared to the isolated case for average density ( 3rd panel) and 
temperature ( bottom panel) are also shown. 


the outputs from 3D simulations. 

5.1. Integrated Intensity Maps 

Using the o utpu ts from our PDR modeling (method 
described in ^3.5|) , we create CO integrated intensity 


maps from the simulation outputs. Note that the local 
emissivity of CO lines does take into account the optical 
depth of an associated PDR layer, of given total col¬ 
umn density that depends on local volume density. To 
illustrate the method, we perform post-processing on the 


colliding case (u re i = lOkm/s) with an impact parameter 
of b = O.bRi and Reoriented mixed fields, our fiducial 
model. 


The diagnostics portions of Figures [18] and [19] show in¬ 
tegrated intensity maps of common observational bands 
of 13 CO(J=2-l and 3-2) and 12 CO(J=8-7) as well as 
the 12 CO(J=8-7)/ 13 CO(J=2-l) line ratio. These maps 
assume a depth of 1 pc in the z direction and a cloud 
distance of 3 kpc. Note also that the adopted abundance 
ratio of 13 C to 12 C is 1/60. We see that the CO emission 
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Fig. 17.— A resolution study comparing time evolution of av¬ 
erage clump density (Top panel) and temperature ( 2nd panel) for 
the fiducial case (see run 3.D.0). Models at the standard resolu¬ 
tion (1024 2 ) are compared with those run with one and two ad¬ 
ditional levels of AMR. The ratios for average density (3rd panel) 
and temperature (bottom panel) compared to the isolated case at 
the particular resolution, are also shown. 


lines trace molecular gas in general, with the higher- J 
lines indeed probing more strongly shocked regions. As 
J increases, higher temperature material is traced, with 
shock fronts of varying strengths being followed. This 
occurs even for low values of nn- While these line emis- 
sivities are most strongly affected by temperature, they 
are also tracers of high density due to the higher critical 
density of the high-J transitions and the dependence of 
nco on tir. Thus, lower temperature, high nn gas is also 
revealed. 

13 CO(J=2-l) and 13 CO(J=3-2) intensity maps show 
fairly similar structures, primarily tracing high-density 
material as well as higher temperature regions. The 
12 CO(J=8-7) map, however, accentuates more strongly 


shocked regions, closely tracing the high-temperature 
dense regions. 

Strongly shocked, high temperature, high density gas 
- potentially a signature of cloud-cloud collisions - pro¬ 
duces the strongest intensity of higher-level lines. Emis- 
sivities at certain J levels as well as their ratios can act 
as diagnostics of a wide range of conditions and poten¬ 
tially determine shock properties and physical conditions 
in the affected gas. 

The final line ratio map further traces high- 
temperature, high-density material, and de-emphasizes 
low-temperature, high-density material. The 12 CO( J=8- 
7)/ 13 CO(J=2-l) line ratio could be an efficient tracer of 
cloud collisions. 

Figure [2T| explores this potential cloud-collision signa¬ 
ture. The average 12 CO(J=8-7)/ 13 CO(J=2-l) line in¬ 
tensity ratio within the clump is calculated and followed 
over time for a set of isolated and colliding cases. From 
these results, we see that this parameter is an excellent 
tracer of cloud collisions. While the clump in the iso¬ 
lated case (once it settles into a relatively stable state) 
retains a value of this intensity ratio of ^ 1 — 10, a clump 
experiencing a GMC collision sees much larger values of 
the line ratio, even reaching > 10 3 for ^ re i = 40 km s -1 . 
Collision velocities as low as v Te \ = 5 km s -1 show an 
excess of a factor of ~ 10 with respect to the isolated 
case. 


5.2. Spectra 


From the simulations, synthetic spectra were created in 
order to provide a more direct comparison with observed 
cloud kinematics. While the initial conditions and 2D 
geometry are fairly idealized, we expect these diagnostic 
methods to be of general use, e.g., once outputs from 
3D simulations are available. Emission line spectra of 
various observational volumes within the simulation box 
for the isolated and colliding fiducial case are shown in 
Figure [22] 

The isolated case shows a narrow velocity range of 
dense gas tracers There is also a relatively strong peak 
in 12 CO(J=8-7) as the cloud pinches in on itself due to 
the presence of in-plane magnetic fields, but these shocks 
occur at low velocities. The chosen volume has little ef¬ 
fect on the spectra as the main features are localized 
around the clump region. The same features are present 
in both lines of sight, with greater asymmetry in the x- 
velocity simply due to the off-center initial position of 
the clump relative to GMC 1. The spectra for the col¬ 
liding case show a much wider velocity spread in each 
of the CO emission lines. In the 20 pc x 20 pc box, the 
emission peaks in 12 CO(J=8-7) at multiple narrow ve¬ 
locity bands correspond to the strongest shocks as seen 
in Fig. 19 The high emissivity feature in 12 CO(J=8-7) 
for the 7/-velocity indicates strong shocks traveling north¬ 
ward around the collision region. The strong features 
present in the 8 pc by 5 pc region but not the clump ma¬ 
terial represent shocks compressing, but not yet propa¬ 
gating through, the clump. These shocks, directed in the 
negative-x and y directions, indicate the collision with 
GMC 2. 


Next we measure line-of-sight velocities and velocity 
gradients in the 8 pc by 5 pc rectangular region around 
the clump in the fiducial simulation at the time of max¬ 
imum compression using the 13 CO(J=2-l) spectra. We 
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Fig. 18.— Time evolution of colliding GMCs at 2.05, 4.01, 4.99, 5.96, and 7.92 Myr (1 level AMR version of run 3.D.0 in Table[2]) Here, 
out-of-plane magnetic fields are near critical (B c \ = 65 p,G, B\ = 40 fi G, and Bo = 10 p,G) while an additional in-plane uniform field of 
B x = 10 fiG is applied throughout the simulation. Furthermore, GMC 2 is offset at b = 0.5Ri. (Row 1 ): Maps of wh with magnetic fields 
represented by grey streamlines. (Row 2): Maps of temperature with black vectors representing velocity. (Row 3): 13 CO( J=2-l) integrated 
intensity maps using PDR-based, temperature and density dependent volume emissivities. (Row 4)'- Similarly derived 13 CO(J=3-2) line 
intensity maps. (Row 5): Similarly derived 12 CO(J=8-7) line intensity maps. (Row 6): 12 CO(J=8-7)/ 13 CO(J=2-l) line intensity ratio 
maps. The advective scalar defining the clump is shown by black or white contour lines, representing the scalar value S = 0.25, 0.5, and 
0.75. 
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Fig. 19.— Zoomed in maps of the clump at t = 4.99 Myr, near the time of maximum compression due to the collision. This is the 
x-directed mixed field geometry case with b = 0.5Ri (2 level AMR version of run 3.D.0 in Table [ 2 ]) Here, the total H-field magnitude is 
near critical while an additional in-plane uniform field of B x = 10 p,G is applied throughout the simulation. (left 4 figures) Maps of density 
(nn)> temperature, R-field magnitude, and velocity gradient magnitude are shown. Grey streamlines indicate magnetic field structure while 
velocities are represented by the black vectors. (right 4 figures ) Maps of 13 CO(J=2-l), 13 CO(J=3-2), and 12 CO(J=8-7) intensity, as well 
as a map of 12 CO( J=8-7)/ 13 CO( J=2-l) line intensity ratio are shown. The advective scalar defining the clump is shown by black or white 
contour lines, representing the scalar value S = 0.25, 0.5, and 0.75. 
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Fig. 20. — A map of the ratio of actual simulation temperature to 
the density-based equilibrium temperature at t = 4.99 Myr for the 
2 level AMR fiducial case. The advective scalar at values S = 0.25, 
0.5, and 0.75 defining the clump is shown by black contour lines. 


compare gradients derived from the total mass distribu¬ 
tion and those derived from the intensities of 13 CO( J=2- 
1), 13 CO(J=3-2) and 12 CO(J=8-7) spectra. Figure [23 


shows the mean velocities and derived velocity gradients 
of the clump material along orthogonal lines of sight. 
The mean gradients are 0.97 and -0.81 km s -1 pc for 
spectra measured along lines of sight perpendicular and 
parallel to the collision direction, respectively. Somewhat 
smaller gradients are derived from lower J CO lines, and 
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Fig. 21.— Average 12 CO(J=8-7)/ 13 CO(J=2-l) line intensity 
ratio (top) over time comparing the effects of collision velocity for 
the R m i x field geometry case (see runs 3.C.l.x). Here, out-of-plane 
magnetic fields are near critical (B c \ = 65 fiG, B\ = 40 juG, and 
Ho ps 10 fiG) while an additional in-plane uniform field of B x = 
10 fiG is applied throughout the simulation. Collision velocities 
i7 re i = 5,10,20,40km s -1 are shown, along with the evolution of 
the isolated GMC with clump. The ratios of the colliding cases 
compared to the isolated case (bottom) are also shown. 


larger gradients from higher J lines. 

Velocity gradients have been measured observation- 
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Fig. 22 .— Synthetic spectra for 13 CO(J=2-l), 13 CO(J=3-2), 
and 12 CO(J=8-7) are shown for ( top six: (a) - (f)) the isolated 
case and ( bottom six: (g) - (l)) colliding case, at t = 4.99 Myr. The 
subplots denote emission analyzed from different volumes (x X y 
values listed below, and assuming 1 pc extent in the z direction) 
in the simulation: a 20 pc X 20 pc box centered on the clump, a 
smaller 5pcx8pc (isolated case) or 8pcx5pc (colliding case) region 
containing the clump, and contribution solely from the clump ma¬ 
terial, defined by the scalar value S > 0.5. The left column shows 
spectra derived from v x (i.e., a view along the collision axis), while 
the right column shows spectra derived from v y (i.e., a view per¬ 
pendicular to the collision axis). Velocity bins of 0.25 kms -1 are 
used. 


ally within IRDCs and GMCs. For example, |Ragan| 


et al. (2012) found velocity gradients of 2.4 and 2.1 
km s _1 pc _r within sub-pc regions of IRDCs G5.85-0.23 
and G24.05-0.22 , base d on observations of NH 3 (1, 1). 


Henshaw et al.| (|2014j) found values of 0.08, 0.07, and 
0.30 km s -i pc _r on ^ 2 pc -1 scales and larger local 
gradients (1.5 2.5 km s -1 pc -1 ) on sub-parsec scales in 
IRDC G035.3900.33, based o n centroid velocities of the 
dense gas tracer N2H + (l-0). Hernandez & Tan (2015) 
derived a mean velocity gradient of 0.24 km s -1 pc -i of 
10 IRDC clumps, based on 13 CO(1-0) emission. 

The gradients seen in our simulated clump are some- 
what larger than thos e observed towards IRDCs by|Her-| 


nandez & Tan (2015), which may indicate these IRDCs 
are not being disturbed kinematically by the kind of col¬ 
lision modeled in our fiducial simulation. However, the 
results of a range of 3D simulations and a wider vari- 




Fig. 23.— Mean velocities of an 8 pc by 5 pc rectangular re¬ 
gion around the clump at t — 4.99 Myr, along the (top) ^-direction 
and ( bottom ) y-direction. Black lines represent line of sight veloc¬ 
ities weighted by the mass distribution of the region. Blue crosses 
are intensity-weighted mean velocities derived from 13 CO(J=2-l) 
spectra of 1 pc wide strips that evenly divide the region along 
the line of sight. Green and red crosses denote similarly calculated 
mean velocities from 13 CO(J=3-2) and 12 CO(J=8-7), respectively. 
The colored dotted lines show the best linear fit to each corre¬ 
sponding set of points, with the value of this gradient displayed in 
parenthesis in the legend. 


ety of viewing angles are needed before more definitive 
conclusions can be drawn from such comparisons. 

6. DISCUSSION AND CONCLUSIONS 

We have explored a wide range of parameter space of 
magnetized GMC-GMC collisions. We performed ideal¬ 
ized 2D simulations from the GMC-scale down to ^ 0.1 
parsec scales, allowing us to study in detail the structure 
of GMCs undergoing collisions. In particular, we focused 
on a clump embedded in a GMC, aiming to isolate the 
effects of various parameters on the evolution of clump 
material. 

We began by developing new heating and cooling func¬ 
tions that depend on density, temperature and extinc¬ 
tion, based on the method of VLBT2013. We combined 
the results from the PDR codes PyPDR and Cloudy 
to create arrays of heating and cooling rates that span 
the atomic to molecular transition, allowing treatment of 
a multi-phase ISM and including the thermal instability 
of warm and cold atomic media. Our heating and cool- 














































































Magnetized GMC Collisions 


21 


ing functions return self-consistent rates from densities 
ranging from nn = 10 -3 — 10 6 cm -3 and temperatures 
ranging from T = 2.7 — 10 5 K. This enabled us to study 
non-equilibrium temperature conditions that are present 
in the shocked material of colliding clouds. Further, we 
similarly derived emissivity arrays for various common 
observational bands of CO, allowing us to simulate syn¬ 
thetic observations via post-processing. 

In terms of MHD simulations, our models tracked an 
initially magnetically subcritical clump embedded within 
a GMC. We investigated different collision velocities, im¬ 
pact parameters, magnetic field strengths and orienta¬ 
tions and their effects on colliding versus isolated GMCs. 
For the maximally supportive out-of-plane -field cases, 
we reported GMC collisions at typical velocities caus¬ 
ing density increases of a factor of ^ 2 — 3, with a 
rebound resulting in relatively lower average densisties 
post-collision. During the collision, average clump tem¬ 
peratures were increased by a factor of up to ^ 10 — 20 
due to shocks dominating the clump material before set¬ 
tling back to ^ 15 — 30 K. Collisions with impact param¬ 
eter between the GMCs produced similar levels of con¬ 
traction with less exaggerated effects for higher impact 
parameters. However, these types of collisions involve 
strongly shearing velocity fields that produced asymme¬ 
tries and more filamentary structure, as well as imparting 
angular momentum to the resulting cloud. 

Mixed-field geometries resulted in relative increases of 
density and temperature at levels similar to the out-of¬ 
plane case. However, late-time behavior of the clouds 
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showed eventual contraction, as material is able to flow 
along the field lines and slowly accumulate onto the 
clump. 

Analysis of CO line emissivities provided a way to track 
shocks of various strengths. In particular, the average 
value of the 12 CO(J=8-7)/ 13 CO(J=2-l) ratio within a 
clump that was undergoing a collision versus one in an 
isolated cloud resulted in differences of a factor of up 
to ~ 10 4 for typical collision velocities. Even slow col¬ 
lisions of = 5 km s -1 showed excesses of a factor of 
^ 10 in this parameter. This may be a useful diagnostic 
signature of cloud collisions. Spectra and velocity gradi¬ 
ents of molecular line emission around dense gas clumps 
may also provide tests of cloud collisions as a triggering 
mechanism for their formation. 

One caveat of the presented models is that all the 
shocks are in the context of ideal MHD. Also, the effects 
of initial GMC turbulence, ambipolar diffusion, star for¬ 
mation, and stellar feedback have not been addressed in 
this study, but are planned in future 3D models. How¬ 
ever, star formation and stellar feedback are not likely to 
be too important in comparing simulation outputs with 
some ISM clouds, such as Infrared Dark Clouds. A more 
complete study of observational diagnostics with compar¬ 
ison to cases in the Galaxy is planned in a subsequent 
paper that will analyze 3D simulations and include initial 
GMC turbulence. 

We thank Fumitaka Nakamura and Shuo Kong for use¬ 
ful discussions. J.C.T. acknowledges NASA Astrophysics 
Theory and Fundamental Physics grant ATP09-0094. 


Jonkheid, B., Faas, F. G. A., van Zadelhoff, G.-J., & van 
Dishoeck, E. F. 2004, A&A, 428, 511 
Kennicutt, Jr., R. C. 1998, ApJ, 498, 541 
Klein, R. I., & Woods, D. T. 1998, ApJ, 497, 777 
Koda, J., Sawada, T., Hasegawa, T., & Scoville, N. Z. 2006, ApJ, 
638, 191 

Kortgen, B., & Banerjee, R. 2015, MNRAS, 451, 3340 
Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304 
Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 
Li, H.-B., Goodman, A., Sridharan, T. K., et al. 2014, Protostars 
and Planets VI, 101 

Liszt, H. S., Burton, W. B., & Xiang, D.-L. 1984, A&A, 140, 303 
Mac Low, M.-M., Smith, M. D., Klessen, R. S., &; Burkert, A. 

1998, Ap&SS, 261, 195 
Matzner, C. D. 2007, ApJ, 659, 1394 
McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 
O’Shea, B. W., Bryan, G., Bordner, J., et al. 2004, ArXiv 
e-prints, astro-ph/0403044 

Ostriker, E. C., Gammie, C. F., &; Stone, J. M. 1999, ApJ, 513, 
259 

Ostriker, J. 1964, ApJ, 140, 1056 

Ragan, S. E., Heitsch, F., Bergin, E. A., & Wilner, D. 2012, ApJ, 
746, 174 

Rollig, M., Ossenkopf, V., Jeyakumar, S., Stutzki, J., & 
Sternberg, A. 2006, A&A, 451, 917 
Rollig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187 
Rosolowsky, E., Engargiola, G., Plambeck, R., & Blitz, L. 2003, 
ApJ, 599, 258 

Sarazin, C. L., & White, III, R. E. 1987, ApJ, 320, 32 
Schoeier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & 
Black, J. H. 2005, VizieR Online Data Catalog, 343, 20369 
Scoville, N. Z., Sanders, D. B., & Clemens, D. P. 1986, ApJ, 310, 
L77 

Smith, B., Sigurdsson, S., & Abel, T. 2008, MNRAS, 385, 1443 

Stark, A. A. 1984, ApJ, 281, 624 

Sternberg, A., & Dalgarno, A. 1989, ApJ, 338, 197 


22 


Wu et al. 


Suwannajak, C., Tan, J. C., & Leroy, A. K. 2014, ApJ, 787, 68 
Takahira, K., Tasker, E. J., & Habe, A. 2014, ApJ, 792, 63 
Tan, J. C. 2000, ApJ, 536, 173 
—. 2010, ApJ, 710, L88 

Tan, J. C., Shaske, S. N., & Van Loo, S. 2013, in IAU Symp. 292, 
Molecular Gas, Dust, and Star Formation, ed. T. Wong & J. 
Ott (Cambridge: Cambridge Univ. Press), 19 
Tasker, E. J., & Tan, J. C. 2009, ApJ, 700, 358 
Tielens, A. G. G. M. 2005, The Physics and Chemistry of the 
Interstellar Medium (Cambridge: Cambridge Univ. Press) 
Tomisaka, K. 2014, ApJ, 785, 24 

Troland, T. H., & Crutcher, R. M. 2008, ApJ, 680, 457 
Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, 
L179 


Van Loo, S., Butler, M. J., & Tan, J. C. 2013, ApJ, 764, 36 
van Loo, S., Falle, S. A. E. G., & Hartquist, T. W. 2010, 
MNRAS, 406, 1260 

van Loo, S., Falle, S. A. E. G., Hartquist, T. W., & Moore, 
T. J. T. 2007, A&A, 471, 213 
Wang, P., & Abel, T. 2008, ApJ, 672, 752 
Wolfire, M. G., McKee, C. F., Hollenbach, D., & Tielens, 

A. G. G. M. 2003, ApJ, 587, 278 
Zuckerman, B., & Evans, II, N. J. 1974, ApJ, 192, L149 



